# -*- coding: utf-8 -*-
"""
Created on Mon Jun 17 11:18:26 2019

@author: Administrator
"""

import cv2
import numpy as np
import matplotlib.pyplot as plt
 
lena = cv2.imread("lena.jpg")#读取目标图片lena.jpg

#1. 以Lena为原始图像，通过OpenCV实现平均滤波，高斯滤波及中值滤波，比较滤波结果。下面是第1题的实现。
Blur_result = cv2.blur(lena, (3,3)) #平均滤波，以3*3的核大小为例
GaussianBlur_result = cv2.GaussianBlur(lena,(3,3),0)  #高斯滤波，以3*3的核大小为例   
medianBlur_result = cv2.medianBlur(lena,3)  #中值滤波，以3的窗口为例   

#显示图形
titles = ['Source', 'Blur', 'GaussianBlur', 'medianBlur']  
images = [lena, Blur_result, GaussianBlur_result, medianBlur_result]  
for i in np.arange(4):  
   cv2.imshow(titles[i],images[i])
cv2.waitKey(0)


#2. 以Lena为原始图像，通过OpenCV使用Sobel及Canny算子检测，比较边缘检测结果。

#彩色转灰色
gray = cv2.cvtColor(lena,cv2.COLOR_BGR2GRAY)

sobel = cv2.Sobel(gray,-1,1,0,ksize=3) #使用sobel算子进行边缘检测,ksize是指核的大小,是对x方向进行边缘检测
laplacian = cv2.Laplacian(gray,-1) #使用laplacian算子进行边缘检测
canny = cv2.Canny(gray, 10, 300)

#显示图形
titles = ['Source', 'sobel', 'laplacian', 'Canny']  
images = [gray, sobel, laplacian, canny]  
for i in np.arange(4):  
   cv2.imshow(titles[i],images[i])
cv2.waitKey(0)


#3. 在OpenCV安装目录下找到课程对应演示图片(安装目录\sources\samples\data)，首先计算灰度直方图，进一步使用大津算法
#进行分割，并比较分析分割结果。

img = cv2.imread("Rice.png")  #导入图片，图片放在程序所在目录

gray_img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) #转换为灰度图
plt.hist(gray_img.ravel(), 256, [0, 256]) #计算灰度直方图
plt.show()


#大津算法进行分割
retval,dst = cv2.threshold(gray_img, 0, 255, cv2.THRESH_OTSU) #使用大津算法进行图像二值化
#显示图形
titles = ['Source', 'dst']  
images = [gray_img, dst]  
for i in np.arange(2):  
   cv2.imshow(titles[i],images[i])
cv2.waitKey(0)

 
#4. 使用米粒图像，分割得到各米粒，首先计算各区域(米粒)的面积、长度等信息，进一步计算面积、长度的均值及方差，
#分析落在3sigma范围内米粒的数量。 

#使用第3题中的灰度图像gray_img，使用局部阈值的大津算法进行图像二值化
dst_1 = cv2.adaptiveThreshold(gray_img,255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY,101, 1)


#全局大津算法，效果较差
#res ,dst = cv2.threshold(gray,0 ,255, cv2.THRESH_OTSU)

element = cv2.getStructuringElement(cv2.MORPH_CROSS,(3, 3)) #形态学去噪
dst_1 = cv2.morphologyEx(dst_1,cv2.MORPH_OPEN,element)  #开运算去噪

im2, contours, hierarchy = cv2.findContours(dst_1,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)  #轮廓检测函数
cv2.drawContours(dst_1,contours,-1,(120,0,0),2)  #绘制轮廓

count=0 #米粒总数
ares_list = []  #面积列表
length_list = []  #长度列表
#遍历找到的所有米粒
for cont in contours:
    ares = cv2.contourArea(cont)#计算包围性状的面积
    length = cv2.arcLength(cont,1)#计算包围性状的长度
    if ares<0:   #过滤面积小于10的形状
        continue
    count+=1    #总体计数加1
    ares_list.append(ares)
    length_list.append(length)
    
    print("{}-blob:{}".format(count,ares),end="  ") #打印出每个米粒的面积

    rect = cv2.boundingRect(cont) #提取矩形坐标

    print("x:{} y:{}".format(rect[0],rect[1]))#打印坐标

    cv2.rectangle(img,(rect[0],rect[1]),(rect[0]+rect[2],rect[1]+rect[3]),(0,0,0xff),1)#绘制矩形
    y=10 if rect[1]<10 else rect[1] #防止编号到图片之外

    cv2.putText(img,str(count), (rect[0], y), cv2.FONT_HERSHEY_COMPLEX, 0.4, (0, 255, 0), 1) #在米粒左上角写上编号

cv2.namedWindow("imagshow", 2)   #创建一个窗口
cv2.imshow('imagshow', img)    #显示原始图片
cv2.waitKey()

count_ares=0 #米粒总落在3sigma范围内（长度）米粒的数量
for i in ares_list:
    if i<=(np.mean(ares_list)+3*np.std(ares_list)) and  i>=(np.mean(ares_list)-3*np.std(ares_list)):
         count_ares+=1 

count_length=0 #米粒总落在3sigma范围内（面积）米粒的数量
for i in length_list:
    if i<=(np.mean(length_list)+3*np.std(length_list)) and  i>=(np.mean(length_list)-3*np.std(length_list)):
         count_length+=1 

print("米粒的总数:{}".format(count)) 
print("面积的均值:{}".format(round(np.mean(ares_list),2)))
print("面积的方差:{}".format(round(np.var(ares_list),2)))
print("落在3sigma范围内（面积）米粒的数量:{}".format(count_ares))
print("长度的均指:{}".format(round(np.mean(length_list),2)))
print("长度的方差:{}".format(round(np.var(length_list),2)))
print("落在3sigma范围内（长度）米粒的数量:{}".format(count_length))


#扩展作业： 
#5. 使用棋盘格及自选风景图像，分别使用SIFT、FAST及ORB算子检测角点，并比较分析检测结果。 
#(可选)使用Harris角点检测算子检测棋盘格，并与上述结果比较。 

#先对棋盘格图片进行检测   
#读入棋盘格图片
img = cv2.imread("left04.jpg")
cv2.imshow('Input Image', img)
cv2.waitKey(0)
	
#(1)SIFT检测
sift = cv2.xfeatures2d.SIFT_create()
keypoints = sift.detect(img, None)
	
# 显示
# 必须要先初始化img2
img2 = img.copy()
img2 = cv2.drawKeypoints(img, keypoints, img2, color=(0,255,0))
cv2.imshow('Detected SIFT keypoints', img2)
cv2.waitKey(0)

#(2)FAST检测
# 注意有_create()后缀
fast = cv2.FastFeatureDetector_create()
keypoints = fast.detect(img, None)
	
# 必须要先初始化img2
img2 = img.copy()
img2 = cv2.drawKeypoints(img, keypoints, img2, color=(0,255,0))
cv2.imshow('Detected FAST keypoints', img2)
cv2.waitKey(0)

#(3)ORB检测
orb = cv2.ORB_create()
keypoints = orb.detect(img, None)
	
# 显示
# 必须要先初始化img2
img2 = img.copy()
img2 = cv2.drawKeypoints(img, keypoints, img2, color=(0,255,0))
cv2.imshow('Detected ORB keypoints', img2)
cv2.waitKey(0)

#(4) Harris角点检测
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
gray = np.float32(gray)
dst = cv2.cornerHarris(gray,2,3,0.04)
img[dst>0.01*dst.max()]=[0,0,255]
cv2.imshow('dst',img)
if cv2.waitKey(0) & 0xff == 27:
    cv2.destroyAllWindows()

#先对风景图片进行检测   
#读入风景图片
img = cv2.imread("timg.jpg")
cv2.imshow('Input Image', img)
cv2.waitKey(0)
	
#(1)SIFT检测
sift = cv2.xfeatures2d.SIFT_create()
keypoints = sift.detect(img, None)
	
# 显示
# 必须要先初始化img2
img2 = img.copy()
img2 = cv2.drawKeypoints(img, keypoints, img2, color=(0,255,0))
cv2.imshow('Detected SIFT keypoints', img2)
cv2.waitKey(0)

#(2)FAST检测
# 注意有_create()后缀
fast = cv2.FastFeatureDetector_create()
keypoints = fast.detect(img, None)
	
# 必须要先初始化img2
img2 = img.copy()
img2 = cv2.drawKeypoints(img, keypoints, img2, color=(0,255,0))
cv2.imshow('Detected FAST keypoints', img2)
cv2.waitKey(0)

#(3)ORB检测
orb = cv2.ORB_create()
keypoints = orb.detect(img, None)
	
# 显示
# 必须要先初始化img2
img2 = img.copy()
img2 = cv2.drawKeypoints(img, keypoints, img2, color=(0,255,0))
cv2.imshow('Detected ORB keypoints', img2)
cv2.waitKey(0)